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Abstract 

We study numerically and analytically a stochastic group selection model in which 
a population of asexually reproducing individuals, each of which can be either al- 
truist or non-altruist, is subdivided into M reproductively isolated groups (demes) 
of size N. The cost associated with being altruistic is modelled by assigning the 
fitness 1 — r, with t € [0, 1], to the altruists and the fitness 1 to the non- altruists. 
In the case that the altruistic disadvantage r is not too large, we show that the 
finite M fluctuations are small and practically do not alter the deterministic results 
obtained for M —* oo. However, for large r these fluctuations greatly increase the 
instability of the altruistic demes to mutations. These results may be relevant to 
the dynamics of parasite-host systems and, in particular, to explain the importance 
of mutation in the evolution of parasite virulence. 
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1 Introduction 



Despite the scarcity of empirical facts supporting group selection as a relevant 
evolutive force in nature, the mathematical problems involved in its modeling 
have kept a recurrent theoretical interest on this controversial theory [1,2]. 
Group selection is based on an analogy between individuals (or genes) and re- 
productively isolated subpopulations, termed demes. If the extinction of demes 
occurs at a rate depending on their composition, then such extinctions will fa- 
vor the existence of individuals that increase the probability of survival of the 
deme they belong to. In the case that these individuals are disfavored by the 
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usual selection at the individual level, group selection will oppose individual 
selection, and so it has been advanced as an explanation for the existence of 
altruistic traits in nature. Such a trait is defined as one that is detrimental to 
the fitness of the individual who expresses it, but that confers an advantage 
on the group of which that individual is a member. 

The standard mathematical framework to study group selection was proposed 
by Levins more than 20 years ago [1,2]. The key ingredients are the differ- 
ential survival probability favoring demes with a large number of altruists, 
and the subsequent recolonization of the extinguished demes by the surviv- 
ing ones. In fact, this is practically the only generally accepted mechanism 
to produce group selection in nature (see [3,4] for an alternative proposal). 
However, the mathematical complexity of Levins' formulation, based on a 
nonlinear integral partial differential equation, as well as the need for too 
restrictive assumptions, have motivated the proposal and study of a variety 
of discrete time versions of Levins' model [5-8]. These analyses have concen- 
trated mainly on the deterministic regime, in which the number of demes M is 
infinite, though the deme size N (i.e., the number of individuals in each deme) 
is finite. In the absence of mutation, the finitude of N is crucial to guarantee 
the fixation through random drift of the altruistic trait within some demes. 
The undesirable feature of considering M finite as well, besides obliterating 
the possibility of an analytical solution to the problem, is that the fluctuations 
occurring during the extinction process will ultimately lead to the complete 
extinction of the population. This is a consequence of the sequential procedure 
that considers first the extinction of demes and then the recolonization of the 
extinct demes by the surviving ones. In this paper we modify the sequential 
extinction-recolonization procedure so as to avoid global extinction, allowing 
thus the numerical and analytical study of the effects of a finite population on 
the steady-state of the metapopulation (i.e., the population of demes). More 
pointedly, once a deme is extinct we immediately assign one of the M — 1 sur- 
viving demes to replace it, although the effective replacement of the extinct 
demes will take place only after all demes have passed the extinction stage. 
This replacement or recolonization occurs simultaneously for all demes. 

Our goal is to study the effects of the fluctuations due to the finitude of 
the population on the stability of the altruistic state predicted by Eshel in the 
deterministic regime [6]. The remainder of this paper is organized as follows. In 
section 2 we describe the events that comprise the life cycles of the individuals 
in the metapopulation and present the stochastic dynamics governing the time 
evolution of the metapopulation. A mean-field recursion equation is derived 
and its validity discussed in section 3. The results of the simulations as well as 
those of the mean-field approximation are presented and analyzed in section 
4. Finally, in section 5 we present some concluding remarks and, in particular, 
point out the relevance of our results to the dynamics of parasite-host systems. 
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2 Model 



The metapopulation is composed of M demes, each of which being composed 
of N haploid, asexually reproducing individuals. An individual can be either 
altruist or non- altruist. The cost associated with being altruistic is modelled 
by assigning the reproductive rate 1 — r, with r G [0, 1], to the altruists and the 
reproductive rate 1 to the non-altruists. The demes are classified according to 
the number of altruistic individuals they have, so that there are N + l different 
types of demes, labeled by the integers % — 0, 1, . . . , N. In each generation the 
metapopulation is described by the vector n = (no,ni, . . . ,un), where n« is 
the number of demes of type i, so that J2i n i — M. The life cycle (i.e., one 
generation) consists of the following events, which will be discussed in detail 
in the sequel: extinction, recolonization, reproduction, and mutation. 



2. 1 Extinction and Recolonization 



Within the differential extinction framework, we define the probability that a 
deme of type i survives extinction, a iy by 



an 



Ul + i/i c ) 



if i < i c 
otherwise, 



(1) 



where i c = 0,1, ... ,N is a parameter measuring the intensity of the group 
selection pressure. The larger the number of altruists in a deme, the larger 
its chance of surviving extinction. Once a deme is extinct, a randomly chosen 
deme among the M — l surviving ones will immediately be assigned to replaced 
it. This contrasts with the standard modelling in which recolonization takes 
places only after all demes have passed the extinction procedure. Hence, given 
n, the probability that a deme of type j changes to a deme of type i, denoted 
by Eij, is simply 

a i + (l-a i )(n i -l)/(M-l) if i = j 
(l-a^m/iM-l) if / / ./• 



= < 



As expected, — 1 Vj. We note that the transition matrix E depends 

on n and so it changes as the population evolves. The conjunction of the 
extinction and recolonization procedures is termed interdemic selection since 
the correlation between the elements belonging to a same column of E yields 
an effective, indirect interaction between the demes. 
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2. 2 Reproduction 



The reproduction process occurs inside the denies and hence is termed intra- 
demic selection. Since the size of the denies is fixed and finite (JV), random 
drift occurs. Following Wright's classical model [9] we assume that the number 
of offspring that an individual contributes to the new generation is propor- 
tional to its relative reproductive rate. Thus, the probability that a deme of 
type j changes to a deme of type i is written as 

w) (1 - wjf' 1 , (3) 




where 



J N-jt w 



is the relative reproductive rate of the subpopulation of altruists in a deme 
of type j. We note that J2iRij — 1 Vj and J2iiRij — Nwj. In the absence of 
mutations, the random drift inherent to the reproduction process will prevent 
the existence of mixed demes, i.e., a deme will be either of type N (only altru- 
ists) or of type (only non-altruists). As a result, a choice of the parameter % c 
different from or N in the definition of the survival probability ctj will have 
practically no effect on the metapopulation evolution. 



2.3 Mutation 



To include mutation into the model, we must descend to the level of the genes 
that determine the characteristics of the individuals. In particular, we assume 
that two alleles, say A or B, at a single locus determine whether a given 
individual is altruist or non- altruist, respectively. Since the replication of a 
gene may not be perfect, we introduce the mutation rate u G [0, 1/2], which 
gives the probability that the allele A mutates to B and vice-versa. Hence the 
probability that a deme of type j changes to a deme of type i due to mutations 
of its members is given by 




u 



i+j-2l 



[1-U) 



N-i-j+2l 



(5) 
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where k = max (0, i + j — N) and l u = min Clearly, J2i Uij = 1 Vj and 
E i iU ij = Nu + j(l-2u). 



2.4 Stochastic dynamics 

Given a population characterized by the vector n which will change into a 
new population characterized by n' due to a generic transition matrix T 
(J2iT{j = 1 Vj), the stochastic dynamics is defined by the conditional prob- 
ability distribution P? (n']n). To evaluate this quantity, it is more convenient 
to introduce the set of integers {bij}, where b^ stands for the number of 
demes of type j that have changed to a deme of type %. Hence rij = J2i bij 
and n\ = so that given the set {%}, the vector n' can be read- 

ily determined. In fact, given rij the conditional probability distribution of 
bj = (b j, b^, ... , bjyj) is a multinomial 

Pt (b,- |n,-) = - f , T t ■ ■ ■ T n7 (6) 

j \ 0\j\ . . . N j\ 

for j = 0, . . . , N. Clearly, the random variables bkj and b h are statistically 
independent for i 7^ j, regardless of the values assumed by the indices k 
and /. We must emphasize that since the transition matrix E, which governs 
the extinction and recolonization procedures, depends explicitly on n, this 
dynamics must be applied in parallel (simultaneously) to all demes. 

The dynamics proceeds as follows. Given the population vector in generation 
t, denoted by n', first we consider the extinction-reco Ionization event and 
generate the conditional probability distributions Pe (bj|n*^ Vj. The choice 
of (N + l) 2 uniformly distributed random numbers allows the determination 
of the set {b^} and, consequently, of the new population vector n'. Next, given 
n' we repeat this procedure for the reproduction event, then generating the 
population vector n". Finally, the same procedure is repeated again for the 
mutation event, leading from n" to n* +1 , which thus completes the life cycle. 



3 Expectations 

The (conditional) expected value of the number of demes of type i after a life 
cycle given the population vector n* in generation t is defined by 

(n* +1 ) = E k i p u (k|l) Pr (l|m) P E (m|n*) (7) 

k,l,m 
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where the conditional probabilities are given by Eq. (6) with the generic tran- 
sition matrix replaced by the specific matrices U, R and E as indicated. Here 
we have used the notation 



N 



3=0 



(8) 



with ki = J2j hj Vi. Moreover, using 



E*i^r(k|l)=E% 



(9) 



we find 



jkl 



(10) 



which, by making explicit the dependence of E on n 1 , can be rewritten as 



(nl +1 } = J2 U i3 R 3k\n t k a k + 



jk 



M — 1 



2(l-Q!,)nf-(l-a fc ) 
. i 



At this stage we can readily derive a mean-field recursion equation for the 
average number of demes of type i. In fact, assuming that the covariance 



Cov (nj,n5) = (njn5>-(n*>(n5> 



'.12) 



vanishes at any t for all pairs and setting (nj) = v\ yield 



^ 1 = Y, u ii R iA v i a *+ 



jk 



M — 1 



. i 



■ (13) 



Thus, rather than studying the evolution of a specific population, in this ap- 
proximation scheme we focus on the evolution of an average population whose 
deme frequencies at each generation are regarded as the average of the deme 
frequencies of an infinite number of populations at that generation. Of course, 
for finite M the covariance can never vanish for all pairs since the ran- 
dom variables n\ are not statistically independent (for instance, they obey the 
normalization condition Y^i n \ — M). However, depending on the values of 
the control parameters r, u and M, either the covariance Cov(rik,ni) or its 
coefficient J2j UijRjk(l — ai) may be sufficiently small so as to validate the 
mean-field equation as a good approximation. 
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The (conditional) covariance after a life cycle given n* in generation t is simply 
given by 

Cov (n*, n$) =-E S ik S jk n\ i ± j (14) 
k 

where SV,- = (URE)ij is the matrix element of the product of the three 
transition matrices. The case i = j corresponds to the variance, Var (nf) = 
Cov (n*, ra*), and yields 

Vax(^) =£S iJfe [l (15) 

We note that these quantities can be readily evaluated within the mean-field 
approach by replacing n\ by its average, v\. Of course, the estimate of the 
magnitude of the fluctuations of the random variables (i — 0, . . . , N) around 
their means is crucial to assess the relevance of the finite M effects. 



4 Analysis of the results 

The quantity of interest is the fraction of altruistic individuals in the metapop- 
ulation in the stationary regime, defined as 

1 N 

P = ttE^ (16) 
iV i=0 

where Yi = rii/M is the frequency of demes of type %. Clearly, X^Y^ = I. 
Interestingly, in the case i c = N, there is a simple relation between p and the 
mean fitness of the metapopulation which is defined as a = ^ZiO-iYi^ namely, 

a = \{l+p). (17) 
To measure the dispersion of the random variable p we introduce the variance 

-l = ^h 2 ^-p 2 as) 

i=0 

which vanishes in the case of an homogeneous metapopulation (Yk = 1 and 
Yi = for % 7^ k), and reaches the maximum value 1/4 in the case that 
the demes are segregated in equal proportions into the two opposite classes 
(Yo = Y N = 1/2). Since p and a 2 v are random variables, we will focus on 
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their average values, denoted by (p) and respectively. In all simulations 
discussed in this work, the symbols represent the averages over 2 10 3 indepen- 
dent experiments. The error bars are calculated by measuring the standard 
deviation of the average results obtained in 50 sets of experiments, each one 
involving 40 independent runs. Moreover, in each run the population is left to 
evolve for 2 10 3 generations and we average over the quantities under analysis 
in the last 100 generations. No significant differences were found for longer 
runs. Throughout our analysis we set i c = N = 10, so that Eq. (17) holds 
true. 

In Figs. 1(a) and 1(b) we present (p) and (up 1 ), respectively, as functions of 
the mutation rate u for r = 0.9 and two representative values for the number 
of demes, M = 10 and M = 100. For u = the altruistic trait always takes 
over the metapopulation, provided that there is at least one altruistic deme 
in the initial state (n° N > 1). Besides the stable fixed point presented in these 
figures, for large r the mean-field recursion equations possess an unstable one, 
p = u, which can be reached by starting the iteration with non-altruistic demes 
only (uq = M). Thus the effect of finite M, as shown by the results of the 
simulations, is to increase the instability of the altruistic regime against mu- 
tations by stabilizing the mean-field unstable fixed point. (We note that the 
mean-field results actually shown the opposite tendency, which indicates the 
failure of the approximation in this matter.) This is expected since in a smaller 
metapopulation, chance plays a greater role, and so deleterious mutations ac- 
cumulate with a higher probability, causing a more rapid decrease in the mean 
fitness of the metapopulation. In the case that the size of the metapopulation 
is not fixed but depends on its mean fitness, this positive feedback, termed 
mutational meltdown, leads rapidly to the extinction of the metapopulation 
[10,11]. The agreement between the mean-field predictions and the simula- 
tions are very good for M = 100, except in the region just after the variance 
maximum. Up to this maximum the population is composed almost exclu- 
sively of altruistic and non-altruistic demes, while beyond it the number of 
altruistic demes decreases very rapidly, the sole source of altruistic individuals 
being the mutations within the non-altruistic demes. Clearly, in this scenario 
we have p = u in agreement with the simulation results. The occurrence of a 
pronounced maximum in indicates the existence of a phenomenon similar 
to the error threshold transition of Eigen's quasispecies model for molecular 
evolution [12]. (The formal similarity between group selection and molecular 
evolution models has already been pointed out in ref. [8].) We note that even 
for M = 10, the mean-field approximation yields very good results for small u. 
Of course, the agreement between theory and simulations is more problematic 
for small M, since in this case the probability that the altruistic demes are 
lost from the metapopulation due solely to fluctuations becomes significant, 
leading to the so-called stochastic escape phenomenon [13,14]. This loss is 
practically irreversible as the altruistic selective disadvantage is too high to 
allow for the production of new altruistic demes. 
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The same quantities, (p) and (dp), are presented in Figs. 2(a) and 2(b), ex- 
cept that the altruistic disadvantage is reduced to r = 0.2. In this case, the 
effects of the finite M fluctuations are almost suppressed as illustrated by the 
agreement between the mean-field and the simulation results. The stochastic 
escape phenomenon is not important in this case since new altruistic demes 
can readily be generated due to the small reproductive disadvantage of the 
altruistic individuals. The results for M — > oo are practically indistinguish- 
able from those obtained in the mean-field approximation for M = 100. We 
have verified that changing the values of i c and N alters the frequency of the 
altruistic gene smoothly, leaving its qualitative dependence on the mutation 
rate unaffected. 

A more direct measure of the finite M fluctuations is presented in Figs. 3 
and 4, which show the variances of the fraction of non-altruistic and altruistic 
demes, Var (Y ) and Var(Yjv), respectively, as functions of the mutation rate 
for r = 0.9. We note that although the mean-field approximation describes 
very well the fluctuations outside the region where the transition between 
the altruistic (p ~ 1) and non- altruistic (p ~ u) regimes takes place, it fails 
badly in that region. This failure seems more pronounced in Fig. 4 because 
the range of u coincides with the transition region, but a similar discrepancy 
occurs in Fig. 3 also, where the variance peaks for small u are completely 
overlooked by the mean-field approximation. The situation for r = 0.2 is 
illustrated in Fig. 5 where we present Var (Yjv) as a function of u. In this case 
the mean-field approximation reproduces very well the behavior pattern of the 
simulation results, except for the heights of the peaks which, as expected, are 
underestimated. The results for Var (Y ) are very similar to those shown in Fig. 
5. Of course, the disagreement between simulation and theory is expected since 
we are trying to estimate the size of the fluctuations using an approximation 
scheme that neglects those very same fluctuations. However, the surprisingly 
good agreement shown in Fig. 3 for u outside the transition region suggests 
that a self-consistent iterative scheme, where the covariances calculated in 
the mean-field approximation are used to improve that approximation, may 
describe successfully the finite M fluctuations. As expected, these variances 
tend to zero as the number of demes increases. 



5 Conclusion 

In this paper we have modified the standard implementation of the group 
selection mechanism, which considers first the extinction of the demes and 
then the recolonization of the extinct demes by the surviving ones [1,5,7,8], by 
assigning a recolonizing deme to each extinct deme immediately after its ex- 
tinction, according to Eq. (2). The actual replacement of the extinct demes is 
carried out simultaneously for all demes following the stochastic prescription 
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given in Eq. (6). This modified extinction-recolonization procedure avoids the 
otherwise inevitable global extinction of the population. We have verified, how- 
ever, that this procedure yields qualitatively similar results to those obtained 
with the standard extinction-recolonization procedure in the case that the 
metapopulation survives global extinction long enough to reach a metastable 
state. 

It is important to mention that, in contrast to its original and very criticized 
ecological motivation [15,16], some concepts borrowed from group selection 
have been successfully applied to describe the evolution of parasite-host sys- 
tems [17,18]. In this case the hosts are associated with the demes, while the 
parasites correspond to the individuals inhabiting the demes. The role of the 
altruistic individuals is played by the less virulent parasites which, by having 
a lower growth rate, increase the survival probability of the host. Migration of 
individuals between demes corresponds to horizontal transmission of parasites. 
The transmission of the parasite between parent and offspring generations is 
termed vertical transmission. Interestingly, a well-known result is that, in a 
population of asexual hosts, parasites with vertical transmission alone cannot 
persist if the infected hosts suffer any fitness cost [19,20]. (This result is read- 
ily recognized as Eshel's [6], although no reference to that author is made in 
the specialized literature of parasite-host systems.) Our finding that at cer- 
tain ranges of the mutation rate (around 0.04 in Fig. 1), virulent parasites 
with vertical transmission alone almost take over the population yields evi- 
dences of the major role played by mutations in the evolution of virulence [21]. 
This dominance becomes more pronounced as the host population decreases. 
A more thorough formulation of parasite-host dynamics through the classi- 
cal, discrete time population genetics formalism used to study group selection 
models is still lacking. Such formulation will certainly help to uncover many 
more similarities, as well as overlapping results, between these two fascinating 
research fields. 
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Figure captions 



Fig. 1(a) Average frequency of altruists as function of the mutation rate for 
r = 0.9, M =10 (A), and M = 100 (O)- The solid and dashed curves are the 
mean-field results for M = 10 and M = 100, respectively. The straight line is 
(p) = u. The error bar is omitted when it is smaller than the symbol size. The 
parameters are i c = N = 10. 

Fig. 1(b) Average variance of the frequency of altruists as function of the 
mutation rate. The parameters and convention are the same as for Fig. 1(a). 

Fig. 2(a) Same as Fig. 1(a) but for r = 0.2. 

Fig. 2(b) Same as Fig. 1(b) but for r = 0.2. 

Fig. 3 Variance of the fraction of non-altruistic demes, Var (Y ), as function 
of the mutation rate for r = 0.9. The convention is the same as for Fig. 1(a). 

Fig. 4 Variance of the fraction of altruistic demes, Var(Y N ), as function of 
the mutation rate for r = 0.9. The convention is the same as for Fig. 1(a). 

Fig. 5 Same as Fig. 4, but for r = 0.2. 
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